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Spectral  Shock  capturing  Techniques 


The  effort  to  use  spectral  methods  to  simulate  flows  with  shock  waves 
is  summarized  in  four  published  papers.  In  ({2])  the  authors  study  uniform 
high  order  spectral  methods  to  solve  multi-dimensional  Euler  equations  for 
gas  dynamics.  Uniform  high  order  spectral  approximations  with  spectral  ac¬ 
curacy  in  smooth  regions  of  solutions  are  constructed  by  introducing  the  idea 
of  the  Essentially  Non-Oscillatory  (ENO)  polynomial  interpolations  into  the 
spectral  methods.  Based  on  the  new  approximations,  nonoscillatory  spectral 
methods  which  possess  the  properties  of  both  upwinding  difference  schemes 
and  spectral  methods  were  proposed.  Numerical  results  are  presented  for  the 
inviscid  Burger’s  equation,  and  for  one  dimensional  Euler  equations  includ¬ 
ing  the  interactions  between  a  shock  wave  and  density  disturbance,  Sod’s 
and  Lax’s  shock  tube  problems,  and  the  blast  wave  problem.  Finally,  the 
interaction  between  a  Mach  3  two  dimensional  shock  wave  and  a  rotating 
vortex  is  simulated. 


In  [3]  we  proved  the  existence  of  one-sided  filters,  for  spectral  Fourier  ap¬ 
proximations  of  discontinuous  functions,  which  can  recover  spectral  accuracy 
up  to  the  discontinuity  from  one  side.  We  also  used  a  least  square  proce¬ 
dure  to  construct  such  a  filter  and  test  it  on  several  discontinuous  functions 
numerically. 


This  paper  was  the  prelude  to  the  full  solution  of  the  Gibbs  phenomenon, 
the  start  of  it  given  in  the  following  two  papers:  It  is  well  known  that  the 
Fourier  series  of  an  analytic  and  periodic  function, truncated  after  2 N  +  1 
terms, converges  exponentially  with  N,  even  in  the  maximum  norm.  It  is  also 
known  that  if  the  function  is  not  periodic,  the  rate  of  convergence  deterio- 
rates;in  particular,  there  is  no  convergence  in  the  maximum  norm, although 
the  functione  is  still  analytic.  This  is  known  as  the  Gibbs  Phenomenon.  In  [4] 
we  show  that  the  first  2N-f  1  Fourier  coefficients  contain  enough  information 
about  the  function,  so  that  an  exponentially  convergent  approximation(  in 
the  maximum  norm)  can  be  constructed.  The  proof  is  a  constructive  one  and 
makes  use  of  the  Gegenbauer  polynomials  C^(z).  It  consists  of  two  steps.  In 
the  first  step  we  show  that  the  first  coefficients  of  the  Gegenbauer  expansion 
(based  on  C^(x)  for  0  <  n  <  m  )  of  any  Z/2  function  can  be  obtained  .within 
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exponential  accuracy,  provided  that  both  X  and  m  are  proportional  to(but 
smaller  than  N).  In  the  second  step  we  construct  the  Gegenbauer  expansion 
based  on  C*(x),  from  the  coefficients  found  in  the  first  step.  We  show  that 
this  series  converges  exponentially  with  N,  provided  that  the  original  function 
is  analytic  (though  periodic).  Thus  we  prove  that  the  Gibbs  phenomenon  can 
be  completely  overcome.  Further  work  had  been  done  in  [5].  In  this  paper 
we  discuss  the  wave-resolution  properties  of  the  Fourier  approximations  of  a 
wave  function  with  discontinuities.  It  is  well  known  that  a  minimum  of  two 
points  per  wave  is  needed  to  resolve  a  periodic  wave  function  using  Fourier  ex¬ 
pansions.  For  Chebyshev  approximations  of  a  wave  function  ,  a  minimum  of 
points  per  wave  is  needed.  Here  we  obtain  an  estimate  for  the  minimum  num¬ 
ber  of  points  per  wave  to  resolve  a  discontinuous  wave  based  on  its  Fourier 
coefficients.  In  our  recent  work  on  overcoming  the  Gibbs  phenomenon,  we 
have  shown  that  the  Fourier  coefficients  of  a  discontinuous  function  contain 
enough  information  to  reconstruct  with  exponential  accuracy  the  coefficients 
of  a  rapidly  converging  Gegenbauer  expansion.  We  therefore  study  the  reso¬ 
lution  properties  of  a  Gegenbauer  expansion  where  both  the  number  of  terms 
and  the  order  increase. 
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Parallel  Psuedospectral  Methods 


We  were  concerned  with  the  use  of  parallel  computers  for  spectral  meth¬ 
ods.  For  that  we  have  checked  Domain  Decomposition  techniques.  In  [C]  We 
have  outlined  a  methodology  for  parallizing  spectral  methods.The  influence 
of  interface  boundary  conditions  on  the  ability  to  parallelize  pseudospectral 
multidomain  algorithms  is  investigated.  Using  the  properties  of  spectral  ex¬ 
pansions,  a  novel  parallel  two-domain  procedure  is  generalized  to  an  arbitrary 
number  of  domains  each  of  which  can  be  solved  on  a  separate  processor. 
This  interface  boundary  condition  considerably  simplifies  influence  matrix 
techniques. 

These  techniques  had  been  improved  in  [7].  A  domain  decomposition 
method  is  examined  to  find  a  time  dependent  parabolic  equation.  The 
method  employs  an  orthogonal  polynomial  collocation  technique  on  multi¬ 
ple  subdomains,  the  subdomain  interfaces  arc  approximated  with  the  aid 
of  a  penalty  method.  The  time  discretization  is  implemented  in  an  ex¬ 
plicit/implicit  finite  difference  method.  The  subdomain  interface  is  approx¬ 
imated  using  an  explicit  Dufort  Frankel  method  while  the  interior  of  each 
subdomain  is  approximated  using  an  implicit  backwards  Euler’s  method. 
The  principal  advantage  to  the  method  is  the  direct  implementation  on  a 
distributed  computing  system  with  a  minimum  of  interprocessor  communi¬ 
cation.  Theoretical  results  are  given  for  Legendre  polynomials  while  compu¬ 
tational  results  are  given  for  Chebyshev  polynomials.  Results  are  given  for 
both  a  single  processor  computer  and  a  distributed  computing  system. 

In  [8]  a  method  is  examined  to  approximate  the  interface  conditions  for 
Chebyshev  polynomial  approximations  to  the  solutions  of  parabolic  prob¬ 
lems,  and  a  smoothing  technique  is  used  to  calculate  the  interface  conditions 
for  a  domain  decomposition  method.  The  method  uses  a  polynomial  of  one 
less  degree  than  the  full  approximation  to  calculate  the  first  derivative  so 
that  interface  values  can  be  calculated  by  using  only  the  adjacent  subdo¬ 
mains.  Theoretical  results  are  given  for  the  consistency  of  the  scheme  and 
practical  results  are  presented.  Computational  results  are  given  for  a  fourth 
order  Runga-Kutta  method  in  two  dimentions  and  compuational  results  are 
given  for  an  explicit/explicit  scheme  in  both  one  and  two  dimentions. 
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Spurious  Frequencies  in  Wake  Flows 


This  is  an  example  that  a  good  computation  can  point  out  problems 
in  wind  tunnel  experiments.  In  [9]  A  detailed  numerical  study  of  two- 
dimensional  flow  past  a  circular  cylinder  at  moderately  low  Reynolds  numbers 
has  been  conducted  using  three  different  numerical  algorithms  for  solving  the 
time-dependent  compressable  Navier-Stokes  equations.  It  was  found  that  if 
the  algorithm  and  associated  boundary  conditions  were  consistent  and  sta¬ 
ble,  then  the  major  features  of  the  unsteady  wake  were  well  predicted.  How¬ 
ever,  it  was  also  found  that  even  stable  and  consistent  boundary  conditions 
could  introduce  additional  periodic  phenomena  riminiscent  of  the  type  seen 
in  privious  wind  -tunnel  experiments.  However,  these  additional  frequen¬ 
cies  were  eliminated  by  formulating  the  boundary  conditions  in  terms  of  the 
characteristic  variables.  An  alalysis  based  on  a  simplified  model  provides  an 
explanation  for  this  behaviour. 

Further  study  of  the  effect  of  boundary  conditions  had  been  carried  out 
in  [1].  The  stability  theory  for  finite  difference  Initial  Boundary- Value  ap¬ 
proximations  to  systems  of  hyperbolic  partial  differential  equations  states 
that  the  exclusion  of  Eigenvalues  and  generalized  eigenvalues  is  a  sufficient 
condition  for  stability.  The  theory,  however,  does  not  discuss  the  nature  of 
numerical  approximations  in  the  presence  of  such  eigenvalues.  In  fact,  as 
was  shown  previously  [9],  for  the  problem  of  vortex  shedding  by  a  2-D  cylin¬ 
der  in  subsonic  flow,  stating  boundary  conditions  in  terms  of  the  primitive 
(non-characteristic)  values  may  lead  to  such  eigenvalues  causing  perturba¬ 
tions  that  decay  slowly  in  space  and  remain  periodic  time.  Characteristic 
formulation  of  the  boundary  conditions  avoided  this  problem.  In  this  paper, 
we  reported  on  a  more  systematic  study  of  the  behavior  of  the  (linearized) 
one-dimensional  gas  dynamic  equations  under  various  sets  of  oscillation  - 
inducing  "legal”  boundary  conditions. 
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Collocation  methods  based  on  Wavelets 


We  were  the  first  to  try  to  construct  differentiation  matrices  by  wavelets. 
In  [10]  Lee  Jameson  discussed  the  notion  of  differentiation  matrices.  The 
numerical  solution  of  a  partial  differential  equation  requires  an  easily  ma¬ 
nipulated  spatial  approximation  to  the  derivative  of  the  unknown  function 
as  well  as  some  method  to  march  forward  in  time.  In  general  one  starts 
from  given  values  of  the  unknown  function,  then  a  finite  dimensional  ap¬ 
proximation,  based  on  those  values,  is  constructed.  This  approximation  is 
differentiated  and  the  result  are  read  at  the  gridpoints.  For  example, in  the 
psuedospectral  Chebyshev  method  for  the  disretization  of  the  equation 

dU(x,t)  _  dF[U(x,t)} 
dt  dx 

One  assumes  that  at  a  given  time  the  values  of  U(xj,t)  are  given  for  some 
points  Xj  =  cos(^),  (j  =  0,N).  Then  one  constructs  the  interpolation 
polynomial  through  those  points  and  differentiate  this  polynomial  to  get  ap¬ 
proximate  values  for  at  the  point  xj.  this  procedure  can  be  viewed 

as  a  transformation  from  N  given  values  (of  the  function)  to  new  N  values 
(approximating  the  derivative.  This  is  the  Chebyshev  Differentiation  Matrix. 
The  numerical  algorithm  therefore  is  simple  and  the  boundary  conditions  can 
be  easily  applied.  It  is  natural  to  ask  whether  one  gain  by  using  wavelets 
instead  of  Chebyshev  polynomials.  .  Since  wavelets  are  well  localized  func¬ 
tions  it  is  reasonable  to  conjecture  that  they  might  represent  steep  gradients 
or  the  development  of  a  shock  with  a  relatively  small  number  of  terms,  con¬ 
sider  a  periodic  function  f(x)  given  on  an  equally  spaced  mesh.  Expand  it  in 
wavelet  expansion  and  use  the  derivative  of  this  expansion  as  an  approxima¬ 
tion  to  the  derivative  of  f(x).  Lee  Jameson,  has  proved  that  approximation 
of  a  periodic  function  /( x)  in  a  wavelet  basis  and  the  differention  of  this 
approximation  yields  nothing  more  than  a  finite  difference  approximation  to 
a  derivative.  The  following  is  an  outline  of  the  proof: 

i)  Given  a  periodic  function  /(x)  let  s  represent  the  periodic  scaling  func¬ 
tion  coefficients  of  this  function  on  the  finest  scale.  This  requires  approx¬ 
imating  the  inner  product  of  f(x)  with  the  scaling  function  on  the  finest 
scale.  The  matrix  representation  of  this  approximation  is  circulent  in  form: 
C  :  f  —*  s,  where  /  represents  /(x)  sampled  on  an  evenly  spaced  grid. 
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ii)  Let  D  be  the  mapping  from  the  scaling  function  coefficients  of  f(x)  to 
the  set  of  scaling  function  coefficients  that  represents  the  derivative  of  f(x): 
D  :  s  — >  s.  Since  f(x)  is  periodic  then  the  matrix  form  of  D  is  circulent  in 
form. 

iii)  All  circulent  martrices  of  the  same  size  commute,  therefore  we  can 
apply  the  operator  D  directly  to  /.  The  operator  D  has  the  effect  of  a  finite 
difference  operator,  and  the  proof  will  be  complete. 

Therefore,  the  effect  of  first  approximating  in  a  wavelet  basis,  then  differ¬ 
entiating  in  this  basis  and  finally  converting  back  from  the  wavelet  basis  to 
the  original  function  is  equal  to  applying  the  appropriate  finite  difference  op¬ 
erator  directly  to  the  equally  spaced  sampled  values  of  the  original  function 
f(x).  The  proof  provides  an  insight  into  the  possibility  of  using  wavelets  for 
solutions  of  PDE’s. 

Wavelets  while  not  more  than  known  finite  difference  schemes 
can  provide  a  mechnism  for  automatic  adaptation  of  the  mesh. 


In  [11]  an  effort  to  use  Galerkin  procedure  had  been  exemined  for  the 
Bergers  equation  using  wavelets  based  on  splines. 
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Study  of  High  Order  Finite  Differences  and  Spectral  Methods. 

Several  issues  concerning  the  applications  of  high  order  methods  were  ex¬ 
amined.  In  [12]  we  outlined  the  direction  of  the  numerical  simulations  of  the 
nect  decade  and  the  typical  problems  associated  with  large  scale  computing. 

In  [13]  we  examined  the  way  to  impose  compatibility  conditions  in  Cheby- 
shev  methods.  Often,  in  solving  an  elliptic  equation  with  Neumann  boundary 
conditions,  a  compatability  condition  has  to  be  imposed  for  well-posedness. 

This  condition  involves  integrals  of  the  forcing  function.  When  pseudospec- 
tral  Chebychev  methods  are  used  to  discretize  the  partial  differential  equa¬ 
tion,  these  integrals  have  to  be  approximated  by  an  appropriate  quadrature 
formula.  The  Gauss-Chebychev  (or  any  variant  of  it  like  the  Gauss- Lobatto) 
formula  cannot  be  used  here  since  the  integrals  under  consideration  do  not 
include  the  weight  function.  A  natural  candidate  to  be  used  in  approximating 
the  integrals  is  the  Clenshaw- Curtis  formula;  however,  we  show  in  this  article 
that  this  is  the  wrong  choice  and  it  may  lead  to  divergence  if  time-dependent 
methods  are  used  to  march  the  solution  to  steady  state.  We  develop, in  this 
paper,  the  correct  quadrature  formula  for  these  problems.  This  formula  takes 
into  account  the  degree  of  the  polynomials  involved.  We  show  that  this  for¬ 
mula  leads  to  a  well  conditioned  Chebychev  approximation  to  the  differential 
equations  and  that  the  compatability  condition  is  automatically  satisfied 

In  [14]  We  study  the  stability  of  spectral  approximations  to  scalar  hy¬ 
perbolic  initial-boundary  value  problems  with  variable  coefficients.  Time  is 
discretized  by  explicit  multi-level  of  Runge-Kutta  methods  of  order  3  (for¬ 
ward  Euler  time  differencing  is  included), and  we  study  spatial  discretizations 
by  spectral  and  pseudospectral  approximations  associated  with  the  general 
family  of  Jacobi  polynomials.  We  prove  that  these  fully  explicit  spectral 
approximations  are  stable  provided  their  time-step,  ,  is  restricted  by  the 
CFL-like  condition,  k  where  N  equals  the  spatial  number  of  degrees  of  free¬ 
dom.  We  give  two  independent  proofs  of  this  result,  depending  on  two  differ¬ 
ent  choices  of  appropriate  -weighted  norms.  In  both  approaches,  the  proofs 
hinge  on  a  certain  inverse  inequality  interesting  for  its  own  sake.  Our  result 
confirms  the  commonly  held  belief  that  the  above  CFL  stability  restriction, 
which  is  extensively  used  in  practical  implementations,  guarantees  the  sta¬ 
bility  (and  hence  the  convergence)  of  fully-explicit  spectral  approximations 
in  the  non-periodic  case. 
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In  [15]  we  discuss  a  new  type  of  schemes:  In  a  previous  paper  we  have  pre¬ 
sented  a  new  method  of  imposing  boundary  conditions  in  the  pseudospectral 
Chebyshev  approximation  of  a  scalar  hyperbolic  equation.  The  novel  idea  of 
the  new  method  is  to  collocate  the  equation  at  the  boundary  points  as  well 
as  in  the  inner  grid  points,  using  the  boundary  conditions  as  penalty  terms. 
In  this  paper  we  extend  the  above  boundary  treatment  to  the  case  of  general 
constant  coefficients  hyperbolic  systems  of  equations,  and  we  provide  error 
estimates  for  the  pseudospectral  Legendre  method.  The  same  scheme  can  be 
implemented  also  in  the  general  (even  nonlinear)  case. 

In  [16]  A  competitive  algorithm,  which  allows  the  computation  of  approx - 
iated  polynomial  solutions  of  advection-diffusion  equations  in  the  square, is 
presented.  The  equation  is  collocated  at  a  special  grid  and  the  correspond¬ 
ing  system  is  solved  by  a  low-cost  preconditioned  iterative  procedure.  The 
method  provides  accurate  results  even  when  the  solution  presents  sharp 
boundary-layers. 

In  [17]  The  authors  examine  psudospectral  matrices. A  standard  way  to 
approximate  the  solution  of  differential  problems  by  algebraic  polynomials, 
is  to  use  collocation  methods  based  on  nodes  related  to  Jacobi  polynomials 
such  as  Chebyshev  or  Legendre  polynomials.  The  matrices  corresponding 
to  the  classical  differential  operators  are  known  to  be  full,  non-symmetric 
and  ill-  conditioned.  We  examine  those  relative  to  the  discretization  of  Neu¬ 
mann  problems  in  one  space  dimension.  In  particular,  we  are  concerned  with 
their  approximation  properties,  their  eigenvalues  and  the  possibility  to  find 
appropriate  preconditioners.  Several  choices  are  available  when  imposing 
boundary  conditions  in  the  approximate  problem.  For  instance,  they  can 
be  either  directly  enforced  or  imposed  in  a  variational  way.  This  results  in 
different  behaviors,  which  can  drastically  affect  the  numerical  effect  of  the 
matrices.  We  analyze  the  different  cases,  pointing  out  the  advantages  and 
drawbacks  in  using  each  strategy. 

A  particular  attention  is  payed  to  the  study  of  preconditioning  matrices. 

In  [18]  a  polynomial  approximation  of  functions  of  matrices  is  considered. 
In  solving  a  mathematical  problem  numerically,  we  frequently  need  to  oper¬ 
ate  on  a  vector  by  an  operator  that  can  be  expressed  as  /(A)  , where  A  is  a 
matrix  [e.g.  exp(A),  A-1]  .  Except  for  very  simple  matrices,  it  is  impractical 
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to  construct  the  matrix  f(A)  explicitly  -  Usually  an  approximation  to  it  is 
used.  This  paper  develops  an  algorithm  based  upon  a  polynomial  approxi¬ 
mation  to  f(A)  .  First  the  problem  is  reduced  to  a  problem  of  approximating 
f(z)  by  a  polynomial  in  z,  where  z  belongs  to  a  domain  D  in  the  complex 
plane  that  includes  all  the  eigenvalues  of  A.  This  approximation  problem 
is  treated  by  interpolating  f(z)  in  a  certain  set  of  points  that  is  known  to 
have  some  maximal  properties.  The  approximation  thus  achieved  is  "almost 
best.”  Implementing  the  algorithm  to  some  practical  problems  is  described. 
Since  a  solutio:.  to  a  linear  system,  an  iterative  solution  algorithm  can  be 
based  on  a  polynomial  approximation  to  /.  We  r’ve  special  attention  to  this 
important  problem. 

In  [19]  We  study  several  algorithms  for  computing  the  Chebyshev  spec¬ 
tral  derivative  and  compare  their  roundoff  error.  For  a  large  number  of 
collocation  points,  the  elements  of  the  Chebyshev  differentiation  matrix,  if 
constructed  in  the  usual  way.are  not  computed  accurately.  A  subtle  cause 
is  found  to  account  for  the  poor  accuracy  when  computing  the  derivative  by 
the  matrix-vector  multiplication  method.  Methods  for  accuracy  computing 
the  elements  of  the  matrix  are  presented  and  we  find  that  if  the  entries  of 
the  matrix  are  computed  accurately,  the  roundoff  error  of  the  matrix-  vec¬ 
tor  multiplication  is  as  small  as  that  of  the  transform-recursion  algorithm. 
Furthermore,  results  of  the  CPU  time  usage  are  shown  for  several  differ¬ 
ent  algorithms  for  computing  the  derivative  by  the  Chebyshev  collocatu  n 
method  for  a  wide  variety  of  two  dimentional  grid  sizes  on  both  an  IBM 
mainframe  and  a  Cray  2  computer.  We  find  that  which  algorithm  is  fastest 
on  a  particular  machine  depends  not  only  on  the  grid  size,  but  also  on  small 
details  of  the  computer  hardware  as  well.  For  most  practical  grid  sizes  used 
in  computation,  the  even-odd  decomposition  algorithm  is  found  to  be  faster 
that  transform-recursion  method. 

In  [20]  The  stability  characteristics  of  various  fourth-  and  sixth  order 
spacial  operators  are  assessed  using  the  theory  of  Gustafsson,Kreiss,  and 
Sundstrom(G-K-S)  for  the  semi  discreet  Initial  Boundary  Value  Problem 
(IBVP).  These  results  are  then  generalized  to  the  fully  discreet  case  using 
a  recently  developed  theory  of  Kriess.  In  all  cases,  favorable  comparisons 
are  obtained  between  G-K-S  theory,  eigenvalue  determination,  and  numer¬ 
ical  simulation.  The  conventional  definition  of  stability  is  then  sharpened 
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to  include  only  those  sp-dial  discretizations  that  are  asymptotically  stable 
(bounded,  Left  Half  Llane  eigenvalues).  It  is  shown  that  many  of  the  higher- 
order  schemes  which  are  G-K-S  stable  arc  not  asymptotically  stable.  A  series 
of  compact  fourth  an  f  sixth-order  schemes,  which  are  both  asymptotically 
and  G-K-S  stable  for  the  scalar  case,  are  then  developed. 
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Hierarchical  optimization 


Prof.  Berger  started  to  develope  an  optimization  component  in  our  pro¬ 
gram.  In  [21]  The  idea  of  hierarchicalgradient  methods  for  optimization  is 
considered.  It  is  shown  that  the  proposed  approach  provides  powerful  means 
to  cope  with  some  global  convergence  problems  characteristic  to  the  classi¬ 
cal  gradient  methods.  Concerning  global  convergence  problems,  four  topics 
are  addressed:  The  detour  "effect”  ,  the  problem  of  multi-scale  problems, 
the  problem  of  highly  ill-conditioned  objective  functions,  and  the  problem  of 
local-minima  traps  related  to  ambiguous  regions  of  attractions. 

The  great  potential  of  hierarchical  gradient  algorithms  is  revealed  through 
a  hierarchical  Gauss-Newton  algorithm  for  unconstrained  nonlinear  least- 
squares  problems.  The  algorithm,  while  maintaining  super-linear  conver¬ 
gence  rate  like  the  common  conjugate  -  gradient  or  quasi-New'ton  methods, 
requires  the  evaluation  of  partial  derivatives  with  respect  to  only  one  vari¬ 
able  on  each  iteration.  This  property  enables  economized  consumption  of 
CPU  time  in  case  the  computer  codes  for  the  derivatives  are  intensive  CPU 
consumers,  e.g.,  when  the  gradient  evaluations  of  ODE  of  PDE  models  are 
produced  by  numerical  differentiation. 

The  hierarchical  Gauss-Newto-  algorithm  is  extended  to  handle  interval 
constraints  on  the  variables  and  its  effectiveness  demonstrated  by  computa¬ 
tional  results. 
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ENO  schemes 


In  [22],  we  have  performed  a  comprehensive  numerical  study  regarding 
the  accuracy  and  convergence  of  ENO  (essentially  non-oscillatory)  schemes 
on  model  problems.  We  proposed  a  modified  ENO  scheme,  with  no  additional 
computational  costs,  which  can  enhance  the  accuracy  in  smooth  regions.  The 
behavior  of  the  modified  ENO  scheme  near  shocks  is  similar  to  the  original 
ENO  schemes. 

In  [23],  we  have  applied  high  order  ENO  schemes  to  two  and  three  di¬ 
mensional  compressible  Euler  and  Navier-Stokes  equations.  Practical  issues, 
such  as  vectorization,  efficiency  of  coding,  cost  comparison  with  other  nu¬ 
merical  methods,  and  accuracy  degeneracy  effects,  are  discussed.  Numerical 
examples  are  provided  which  are  representative  of  computational  problems 
of  current  interest  in  transition  and  turbulence  physics.  These  require  both 
non-oscillatory  shock  capturing  and  high  resolution  for  detailed  structures  in 
the  smooth  regions  and  demonstrate  the  advantage  of  ENO  schemes. 

In  [24]  we  have  proposed  a  new  numerical  approach  to  solve  systems  of 
conservation  laws  of  mixed  hyperbolic-elliptic  type,  by  a  generalization  of  the 
hyperbolic  flux  splitting  methods.  The  possibly  elliptic  flux  is  split  into  a  sum 
of  two  hyperbolic  fluxes  with  only  positive/negative  eigenvalues.  Upwinding 
hyperbolic  solvers  can  thus  be  applied  to  each  of  them  separately.  Equipped 
with  time  splitting  this  procedure  is  stable.  Numerical  results  are  shown 
to  illustrate  the  behavior  of  this  method  on  the  van  der  Waals  equation  in 
gas  dynamics.  We  observe  convergence  with  good  resolution  to  admissible 
solutions  containing  phase  transitions  for  Riemann  problems.  We  also  tested 
more  general  initial  conditions. 

Iu  [25],  we  discussed  about  a  framework  to  use  uniformly  high  order  spec¬ 
tral  method  for  solving  multi-dimensional  Euler  equations  of  compressible 
flow.  Uniform  high  order  approximation  with  spectral  accuracy  in  smooth 
regions  of  the  solution  is  coupled  to  a  non-oscillatory  shock  transition  through 
the  introduction  of  ENO  interpolation  idea  into  the  spectral  methods.  The 
proposed  numerical  procedure  possesses  the  properties  of  both  upwinding  fi¬ 
nite  difference  schemes  and  spectral  methods.  Test  problems  include  Burgers 
equations,  one  and  two  dimensional  Euler  equations  including  shock-vortex 
and  shock-wave  interactions. 
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In  [26],  we  discuss  the  application  of  high  order  compact  finite  difference 
methods  for  shock  calculations.  The  main  idea  is  the  definition  of  a  local 
mean  which  serves  as  a  reference  for  introducing  a  local  nonlinear  limiting 
to  control  spurious  numerical  oscillations  while  keeping  the  formal  accuracy 
of  the  scheme.  For  scalar  conservation  laws,  the  resulting  schemes  can  be 
proven  total  variation  stable  in  one  space  dimension  and  maximum  norm 
stable  in  multi  space  dimensions.  Numerical  examples  are  shown  to  verify 
accuracy  and  stability  of  such  schemes  for  problems  containing  shocks. 

In  [27],  we  applied  the  discontinuous  Galerkin  finite  element  method  (the 
second  order  version)  to  2D  Euler  equations  of  gas  dynamics.  Test  problems 
include  the  standard  shock  reflection,  backward  step  problem,  and  double 
Mach  reflection.  Different  choices  of  limiters  are  explored. 
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